%Calculates optimal savings rate 
function[value]=xstar_int(alpha,rho,t,m,maxmatch,maxcont, mint)

    % Note: m1 is a nuisance parameter that BFP set to 0 (see their fn 24)
    % Define virtual income above the kink

    
    %Max max log(x(1+m)+alpha) +log(1-(1-t)x)+mint*x below the kink
    A=[1;-1];
    B=[maxcont;0];
    options = optimoptions('fmincon','Display','off');

    %xstar_ifbelowM = fmincon(@(x) -1*( rho*log(x*(1+m) + alpha) + log(1-(1-t)*x) + mint*x*(1+m)),0.05,A,B,[],[], [],[],[], options);
    %xstar_ifbelowM = fmincon(@(x) -1*( (x*(1+m) + alpha)^rho *(1-(1-t)*x) + mint*x*(1+m)),0.05,A,B,[],[], [],[],[], options);
    %xstar_ifaboveM = fmincon(@(x) -1*( rho*log(x + maxmatch*m + alpha) + log(1-(1-t)*x) + mint*x*(1+m)),0.05,A,B,[],[], [],[],[], options);
    %xstar_ifaboveM = fmincon(@(x) -1*( (x + maxmatch*m + alpha)^rho * (1-(1-t)*x) + mint*x*(1+m)),0.05,A,B,[],[], [],[],[], options);
    %xstar = xstar_ifbelowM.*(xstar_ifbelowM<=maxmatch.*(1+m)) + xstar_ifaboveM.*(xstar_ifaboveM>maxmatch.*(1+m))...
    %    + maxmatch.*(1+m).*(xstar_ifbelowM>maxmatch.*(1+m)).*(xstar_ifaboveM<=maxmatch.*(1+m));

    xstar=fmincon(@(x) -1*(utility_ev(x, maxmatch, m, alpha, t,rho,maxcont,mint)), 0.05,A,B,[],[], [],[],[], options);
    xstar=max(xstar,0);
  
    % Now, we account for the fact that the personal contribution
    % must be a discrete value (i.e. 1%, 2%, etc.)
    x_top=ceil(xstar.*100)./100;
    x_bot=floor(xstar.*100)./100;
 
    % Can't contribute <0 or  > maxcont (likely redundant)
    x_top=max(x_top,0);
    x_top=min(maxcont,x_top);
    x_bot=max(x_bot,0);
    x_bot=min(maxcont,x_bot);

    
    % The employee chooses the round number that gives the highest utility
    %utility_top = ((rho.*log(x_top*(1+m)+alpha) + log(1-(1-t).*x_top)+mint*x_top*(1+m)).*...
    %    (x_top<=maxmatch)... 
    %    +((rho.*log(x_top+maxmatch*m+alpha)+log(1-(1-t).*x_top)+mint*(x_top+maxmatch*m)).*(x_top>maxmatch)));

    
    
    %utility_bot = ((rho.*log(x_bot*(1+m)+alpha) + log(1-(1-t).*x_bot)+mint*x_bot*(1+m)).*...
    %    (x_bot<=maxmatch))... 
    %    +((rho.*log(x_bot+maxmatch*m+alpha)+log(1-(1-t).*x_bot)+mint*(x_bot+maxmatch*m)).*((x_bot>maxmatch)));
    
    %utility_bot = (((x_bot*(1+m)+alpha).^rho.*(1-(1-t).*x_bot)+mint*x_bot*(1+m)).*...
    %    (x_bot<=maxmatch))... 
    %    +(((x_bot+maxmatch*m+alpha).^rho.*(1-(1-t).*x_bot)+mint*(x_bot+maxmatch*m)).*((x_bot>maxmatch)));
    
    utility_top = utility_ev(x_top,maxmatch, m, alpha, t,rho,maxcont,mint) ;
    utility_bot = utility_ev(x_bot,maxmatch, m, alpha, t,rho,maxcont,mint) ;
    
    %And the winner is:
    value=(xstar>=0).*( (x_top).*(utility_top>=utility_bot) + (x_bot).*(utility_top<utility_bot));
    
    
    
    
    %To make the units of the externality comparable to what we did before,
    %We calculate EV relative to the situation where each person
    %chooses optimally with no internality
    %this function calculates that x*
    
    function[value]=xstar_ev(alpha,rho,t,m,maxmatch,maxcont)
    %Max max log(x(1+m)+alpha) +log(1-(1-t)x)+mint*x below the kink
    A2=[1;-1];
    B2=[maxcont;0];
    options2 = optimoptions('fmincon','Display','off');
    
    xstar2 = fmincon(@(x) -1*( rho*log(x + m.*min(x,maxmatch) + alpha) + log(1-(1-t)*x)),0.05,A2,B2,[],[], [],[],[], options2);
    xstar2=max(xstar2,0);
    
    % Now, we account for the fact that the personal contribution
    % must be a discrete value (i.e. 1%, 2%, etc.)
    x_top2=ceil(xstar2.*100)./100;
    x_bot2=floor(xstar2.*100)./100;
 
    % Can't contribute <0 or  > maxcont (likely redundant)
    x_top2=max(x_top2,0);
    x_top2=min(maxcont,x_top2);
    x_bot2=max(x_bot2,0);
    x_bot2=min(maxcont,x_bot2);

    % The employee chooses the round number that gives the highest utility
    %utility_top = ((rho.*log(x_top*(1+m)+alpha) + log(1-(1-t).*x_top)+mint*x_top*(1+m)).*...
    %    (x_top<=maxmatch)... 
    %    +((rho.*log(x_top+maxmatch*m+alpha)+log(1-(1-t).*x_top)+mint*(x_top+maxmatch*m)).*(x_top>maxmatch)));
    utility_top2 = rho*log(x_top2 + m.*min(x_top2,maxmatch) + alpha) + log(1-(1-t)*x_top2);
    utility_bot2 = rho*log(x_bot2 + m.*min(x_bot2,maxmatch) + alpha) + log(1-(1-t)*x_bot2);
    
    %utility_bot = ((rho.*log(x_bot*(1+m)+alpha) + log(1-(1-t).*x_bot)+mint*x_bot*(1+m)).*...
    %    (x_bot<=maxmatch))... 
    %    +((rho.*log(x_bot+maxmatch*m+alpha)+log(1-(1-t).*x_bot)+mint*(x_bot+maxmatch*m)).*((x_bot>maxmatch)));
   
    %And the winner is:
    value=(xstar2>=0).*( (x_top2).*(utility_top2>=utility_bot2) + (x_bot2).*(utility_top2<utility_bot2));
    
    end


    function [value]=utility_ev(x, maxmatch, m, alpha, t,rho,maxcont,mint)
        
        % Define virtual income above the kink

        x_matched = x + m.*min(x,maxmatch);
       
        % Compute EV. There are four cases to consider here, because x*(theta)
        % and the default can both be above of below the maximum match rate.
        xstar_norm=xstar_ev(alpha,rho,t,m,maxmatch,maxcont);
        xstar_norm_matched=xstar_norm + m.*min(xstar_norm,maxmatch);
        
        value = (((x_matched+alpha)./(xstar_norm_matched+alpha)).^rho).*...
            (1-(1-t).*x)+ ((1-t).*xstar_norm_matched)- 1 +mint.*x_matched ;

    end

end

